{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_id": "00000-fc24b67c-a12f-41d8-8c3e-c9e888bb9777",
    "output_cleared": false
   },
   "source": [
    "# Radio-Pharmaceutical Therapy (RPT) with 177Lu\n",
    "\n",
    "Use the output of the simulation named 'rpt_main1.py'. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:33.600422Z",
     "start_time": "2024-09-16T14:46:33.593753Z"
    },
    "cell_id": "00001-882fa32f-4a7f-4f36-bb4f-22a941cc781b",
    "execution_millis": 185,
    "execution_start": 1605009448424,
    "output_cleared": false,
    "source_hash": "1a4b3bc8"
   },
   "outputs": [],
   "source": [
    "# %matplotlib notebook\n",
    "\n",
    "from exercices_helpers import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:33.913947Z",
     "start_time": "2024-09-16T14:46:33.809837Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CT image size = (124, 72, 190)   spacing = (4.0, 4.0, 4.0)   origin=(-248.48800659179688, -104.93299865722656, -378.625)\n"
     ]
    }
   ],
   "source": [
    "# read the ct image and convert to np array\n",
    "data_folder = Path(\"./data\")\n",
    "ct = sitk.ReadImage(data_folder / \"ct1-4mm.nii.gz\")\n",
    "ct_array = sitk.GetArrayFromImage(ct)\n",
    "print(f'CT image size = {ct.GetSize()}   spacing = {ct.GetSpacing()}   origin={ct.GetOrigin()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:34.058811Z",
     "start_time": "2024-09-16T14:46:34.009984Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dose image size = (128, 128, 209)   spacing = (4.41806, 4.41806, 4.41806)   origin=(-279.3589934082031, -316.61400134277346, -459.922)\n"
     ]
    }
   ],
   "source": [
    "# Consider the dose map file\n",
    "output_folder = Path('./output')\n",
    "dosemap_filename = output_folder / \"ctdose_dose.mhd\"\n",
    "dosemap = sitk.ReadImage(dosemap_filename)\n",
    "print(f'Dose image size = {dosemap.GetSize()}   spacing = {dosemap.GetSpacing()}   origin={dosemap.GetOrigin()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:34.192621Z",
     "start_time": "2024-09-16T14:46:34.171939Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resampled dose image size = (124, 72, 190)   spacing = (4.0, 4.0, 4.0)   origin=(-248.48800659179688, -104.93299865722656, -378.625)\n"
     ]
    }
   ],
   "source": [
    "# resample dose map like ct\n",
    "resampled_dose = resample_image_like(dosemap, ct)\n",
    "print(\n",
    "    f'Resampled dose image size = {resampled_dose.GetSize()}   spacing = {resampled_dose.GetSpacing()}   origin={resampled_dose.GetOrigin()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:34.365006Z",
     "start_time": "2024-09-16T14:46:34.319615Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dose map output/ctdose_dose.mhd : min = 0.0    max = 1.1527506300908177e-06\n",
      "Dose map output/ctdose_dose.mhd : min = 0.0    max = 70677.47867197175 Gy\n",
      "Dose map output/ctdose_dose.mhd : mean = 36.780549875776394 Gy\n"
     ]
    }
   ],
   "source": [
    "# convert as a numpy array\n",
    "dosemap_array = sitk.GetArrayViewFromImage(resampled_dose)\n",
    "print(f'Dose map {dosemap_filename} : min = {np.min(dosemap_array)}    max = {np.max(dosemap_array)}')\n",
    "\n",
    "# scale to the injected activity, according to the activity used to compute the dose map\n",
    "injected_activity = 7.4e9\n",
    "simulated_activity = 1e5\n",
    "s = injected_activity / simulated_activity\n",
    "dosemap_array = dosemap_array * s\n",
    "\n",
    "# integrated activity (considering 177Lu decay only)\n",
    "half_life_sec = 6.647 * 3600 * 24  # about 6.647 days\n",
    "k = np.log(2) / half_life_sec\n",
    "dosemap_array = dosemap_array / k\n",
    "print(f'Dose map {dosemap_filename} : min = {np.min(dosemap_array)}    max = {np.max(dosemap_array)} Gy')\n",
    "print(f'Dose map {dosemap_filename} : mean = {np.mean(dosemap_array)} Gy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:42.997357Z",
     "start_time": "2024-09-16T14:46:42.822717Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwsAAAGyCAYAAACr/FguAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfn0lEQVR4nO3de1zUZd7/8TegA6jNmCmgGyplqXgWE6ezxjoZ22Zat5YpmdnqjW7Cpma5aNZKax5Liu2gtL/N9bB3uZu2KGFqJWqh3Kmp28HCXRvUTCZNQWF+f3jPN0dHZZgZGPD1fDy+j5zv95prrhk+6Xz4XNf1DXE6nU4BAAAAwDlCa3sAAAAAAIITyQIAAAAAj0gWAAAAAHhEsgAAAADAI5IFAAAAAB6RLAAAAADwiGQBAAAAgEckCwAAAAA8IlkAAAAA4BHJAgAAAACPSBYQcBs3btTdd9+tVq1aKSQkRCtXrrzkc9avX6+ePXsqPDxc7dq1U05OTsDHicsD8YhgQjwCCHYkCwi448ePq1u3bsrKyqpS+3379ik5OVl9+/ZVUVGRJkyYoEcffVRr1qwJ8EhxOSAeEUyIRwDBLsTpdDprexC4fISEhOidd97RwIEDL9hm8uTJWr16tXbu3GmcGzp0qI4eParc3NwaGCUuF8QjggnxCCAYNajtAdRXJ0+eVHl5uV/6MplMioiI8EtfdUFBQYGSkpLcztlsNk2YMOGCzykrK1NZWZnxuLKyUkeOHNFVV12lkJCQQA0VdZjr9ySVlZUXbUc8oiYQjwgmTqdTP/74o1q1aqXQ0JqbhMJ3pyDlhN+dOHHCGRMT45TklyMmJsZ54sSJ2n5bfiHJ+c4771y0zXXXXeecOXOm27nVq1c7JTl/+uknj8+ZNm2a3z5vjsvreO2114hHjqA5iEeOYDr2799/0Xj0p9r87tSmTRuPffz3f/+3Mbb//u//djZr1szZuHFj56BBg5x2u92tj2+//dZ51113OSMjI50tWrRwPvHEE85Tp065tfnggw+cPXr0cJpMJue1117rXLx48XljWbhwobNNmzbO8PBwZ+/evZ1btmyp3gfqR1QWAqC8vFx2u1379++X2Wz2qS+Hw6HY2FiVl5eTIV/ElClTlJ6ebjwuLS1V69at/fIzQP3k+n8rMjLS730Tj/AW8Yhg4orHK664osZesza/O33yySeqqKgwHu/cuVO//OUvdf/990uS0tLStHr1aq1YsUIWi0Xjxo3ToEGD9PHHH0uSKioqlJycrJiYGG3atEnfffedRowYoYYNG2rmzJmSfl5vNGbMGL311lvKz8/Xo48+qpYtW8pms0mSli1bpvT0dGVnZysxMVHz58+XzWbT3r17FRUV5dNn4guShQC64oorfP4fzXkZLimJiYlRSUmJ27mSkhKZzeYL/kMaHh6u8PDw886bzWb+McRFXWoaBvGImkQ8IpjUxjQ1s7mRzOZGPvZy2qvWLVq0cHv8/PPP69prr9Vtt92m0tJSvfHGG1qyZIn69esnSVq8eLE6duyozZs3q0+fPlq7dq0+//xzvf/++4qOjlb37t317LPPavLkyZo+fbpMJpOys7MVFxenOXPmSJI6duyojz76SPPmzTOShblz52r06NEaOXKkJCk7O1urV6/WokWL9OSTT/r4mVQfuyEFkNPp9MtxubFarcrPz3c7l5eXJ6vVWksjwuWMeEQwIR5R/53203GmwnD2cfbanQspLy/XX/7yFz3yyCMKCQlRYWGhTp065bZWqEOHDmrdurUKCgoknVlL1KVLF0VHRxttbDabHA6Hdu3aZbTxtN7I1Ud5ebkKCwvd2oSGhiopKcloU1tIFhBwx44dU1FRkYqKiiSdKcUVFRWpuLhY0pkS+YgRI4z2Y8aM0ddff61JkyZpz549evnll7V8+XKlpaXVxvBRz7ji8bPPPpMkffvtt8Qjag3xCARObGysLBaLcWRmZl7yOStXrtTRo0f18MMPS5LsdrtMJpOaNm3q1i46Olp2u91oc3ai4LruunaxNg6HQydOnNDhw4dVUVHhsY2rj9pCshBAVBbO+PTTT9WjRw/16NFDkpSenq4ePXooIyNDkvTdd98Z/zBKUlxcnFavXq28vDx169ZNc+bM0euvv26U6QBfuOLxlltukSQ99dRTxCNqDfEInMt/lYX9+/ertLTUOKZMmXLJV3/jjTc0YMAAtWrVys/vq+5izUIA+ePLfn1IFm6//faLvg9Pdx+9/fbbtX379gCOCpcrVzw6HA5ZLBaVlpa6zdsmHlGTiEfgXD9/2fetD+/X5Xz77bd6//339fbbbxvnYmJiVF5erqNHj7pVF0pKShQTE2O02bp1q1tfrrVFZ7e52HqjsLAwhYWFeWzj6qO2UFkAAADAZW/x4sWKiopScnKycS4hIUENGzZ0Wyu0d+9eFRcXG2uFrFarduzYoYMHDxpt8vLyZDabFR8fb7S52Hojk8mkhIQEtzaVlZXKz8+v9TVJVBYCiMoCAACANyrke2Wh4tJNzlFZWanFixcrJSVFDRr8/PXYYrFo1KhRSk9PV7NmzWQ2mzV+/HhZrVb16dNHktS/f3/Fx8dr+PDhmjVrlux2u6ZOnarU1FRjJ7IxY8Zo4cKFmjRpkh555BGtW7dOy5cv1+rVq43XSk9PV0pKinr16qXevXtr/vz5On78uLE7Um0hWQggkgUAAABv+G8akjfef/99FRcX65FHHjnv2rx58xQaGqrBgwerrKxMNptNL7/8snE9LCxMq1at0tixY2W1WtW4cWOlpKRoxowZRhvXeqO0tDQtWLBAV1999XnrjYYMGaJDhw4pIyNDdrtd3bt3V25u7nmLnmtaiJNvo37nmnt66NAhv9xYpEWLFufNY8XFXWj+L+BSkzFCPOJSiEcEk9qIkZ9f818ym327R5XD8aMsluuJcT+hshBAVBYAAAC8UTuVBVwYyUIAkSwAAAB4g2Qh2LAbEgAAAACPqCwEEJUFAAAAb1SoOrsZnd8H/IVkIYBIFgAAALxRO1un4sKYhgQAAADAIyoLAURlAQAAwBsscA42JAsBRLIAAADgDZKFYMM0JAAAAAAeUVkIICoLAAAA3qCyEGxIFgKIZAEAAMAb7IYUbJiGBAAAAMAjKgsBRGUBAADAG0xDCjYkCwFEsgAAAOANkoVgwzQkAAAAAB5RWQggKgsAAADeoLIQbEgWAohkAQAAwBskC8GGaUgAAAAAPKKyEGBUBgAAAKqK+ywEG5KFAGIaEgAAgDeYhhRsmIYEAAAAwCMqCwFEZQEAAMAbVBaCDclCAJEsAAAAeINkIdgwDQkAAACAR1QWAojKAgAAgDeoLASbOlNZyMzM1A033KArrrhCUVFRGjhwoPbu3evW5uTJk0pNTdVVV12lJk2aaPDgwSopKXFrU1xcrOTkZDVq1EhRUVGaOHGiTp92D6r169erZ8+eCg8PV7t27ZSTk1OtMbuSBV8PAACAy4Nr61RfDrZO9ac6kyxs2LBBqamp2rx5s/Ly8nTq1Cn1799fx48fN9qkpaXp3Xff1YoVK7RhwwYdOHBAgwYNMq5XVFQoOTlZ5eXl2rRpk958803l5OQoIyPDaLNv3z4lJyerb9++Kioq0oQJE/Too49qzZo1Nfp+AQAAgNpWZ6Yh5ebmuj3OyclRVFSUCgsLdeutt6q0tFRvvPGGlixZon79+kmSFi9erI4dO2rz5s3q06eP1q5dq88//1zvv/++oqOj1b17dz377LOaPHmypk+fLpPJpOzsbMXFxWnOnDmSpI4dO+qjjz7SvHnzZLPZvBoz05AAAAC8USHfKwNUFvypzlQWzlVaWipJatasmSSpsLBQp06dUlJSktGmQ4cOat26tQoKCiRJBQUF6tKli6Kjo402NptNDodDu3btMtqc3YerjasPT8rKyuRwONwOiWlIAAAA3vF1CpI/1jzgbHUyWaisrNSECRN00003qXPnzpIku90uk8mkpk2burWNjo6W3W432pydKLiuu65drI3D4dCJEyc8jiczM1MWi8U4YmNjfX6PAAAAQG2rk8lCamqqdu7cqaVLl9b2UCRJU6ZMUWlpqXHs379fEpUFAAAA71BZCDZ1Zs2Cy7hx47Rq1Spt3LhRV199tXE+JiZG5eXlOnr0qFt1oaSkRDExMUabrVu3uvXn2i3p7Dbn7qBUUlIis9msyMhIj2MKDw9XeHj4eedZswAAAOAN125IvvYBf6kzlQWn06lx48bpnXfe0bp16xQXF+d2PSEhQQ0bNlR+fr5xbu/evSouLpbVapUkWa1W7dixQwcPHjTa5OXlyWw2Kz4+3mhzdh+uNq4+AAAAgMtFnakspKamasmSJfr73/+uK664wlhjYLFYFBkZKYvFolGjRik9PV3NmjWT2WzW+PHjZbVa1adPH0lS//79FR8fr+HDh2vWrFmy2+2aOnWqUlNTjcrAmDFjtHDhQk2aNEmPPPKI1q1bp+XLl2v16tVej5nKAgAAgDe4KVuwqTPJwiuvvCJJuv32293OL168WA8//LAkad68eQoNDdXgwYNVVlYmm82ml19+2WgbFhamVatWaezYsbJarWrcuLFSUlI0Y8YMo01cXJxWr16ttLQ0LViwQFdffbVef/11r7dNlUgWAAAAvEOyEGzqTLJQlS/NERERysrKUlZW1gXbtGnTRu+9995F+7n99tu1fft2r8cIAAAA1Cd1Jlmoi6gsAAAAeIPKQrAhWQggkgUAAABvkCwEmzqzGxIAAACAmkVlIYCoLAAAAHiD+ywEG5KFACJZAAAA8MZpSWF+6AP+wjQkAAAAAB5RWQggKgsAAADeoLIQbKgsBJArWfD1AAAAuDyc9tPhnf/85z966KGHdNVVVykyMlJdunTRp59+alx3Op3KyMhQy5YtFRkZqaSkJH3xxRdufRw5ckTDhg2T2WxW06ZNNWrUKB07dsytzWeffaZbbrlFERERio2N1axZs84by4oVK9ShQwdFRESoS5cul7w/WKCRLAAAAOCy9cMPP+imm25Sw4YN9c9//lOff/655syZoyuvvNJoM2vWLL344ovKzs7Wli1b1LhxY9lsNp08edJoM2zYMO3atUt5eXlatWqVNm7cqMcee8y47nA41L9/f7Vp00aFhYV64YUXNH36dL366qtGm02bNumBBx7QqFGjtH37dg0cOFADBw7Uzp07a+bD8IBpSAFGZQAAAKCqan43pD/+8Y+KjY3V4sWLjXNxcXHGn51Op+bPn6+pU6fqnnvukST9+c9/VnR0tFauXKmhQ4dq9+7dys3N1SeffKJevXpJkl566SXdddddmj17tlq1aqW33npL5eXlWrRokUwmkzp16qSioiLNnTvXSCoWLFigO++8UxMnTpQkPfvss8rLy9PChQuVnZ3t06dSXVQWAohpSAAAAN7w3zQkh8PhdpSVlXl8xX/84x/q1auX7r//fkVFRalHjx567bXXjOv79u2T3W5XUlKScc5isSgxMVEFBQWSpIKCAjVt2tRIFCQpKSlJoaGh2rJli9Hm1ltvlclkMtrYbDbt3btXP/zwg9Hm7NdxtXG9Tm0gWQAAAEC9ExsbK4vFYhyZmZke23399dd65ZVXdN1112nNmjUaO3asfvvb3+rNN9+UJNntdklSdHS02/Oio6ONa3a7XVFRUW7XGzRooGbNmrm18dTH2a9xoTau67WBaUgBxG5IAAAA3jgt33+XfaaysH//fpnNZuNseHi4x9aVlZXq1auXZs6cKUnq0aOHdu7cqezsbKWkpPg4lrqPykIAMQ0JAADAG/6bhmQ2m92OCyULLVu2VHx8vNu5jh07qri4WJIUExMjSSopKXFrU1JSYlyLiYnRwYMH3d/J6dM6cuSIWxtPfZz9Ghdq47peG0gWAAAAcNm66aabtHfvXrdz//rXv9SmTRtJZxY7x8TEKD8/37jucDi0ZcsWWa1WSZLVatXRo0dVWFhotFm3bp0qKyuVmJhotNm4caNOnTpltMnLy1P79u2NnZesVqvb67jauF6nNpAsBBCVBQAAAG9U+OmourS0NG3evFkzZ87Ul19+qSVLlujVV19VamqqJCkkJEQTJkzQc889p3/84x/asWOHRowYoVatWmngwIGSzlQi7rzzTo0ePVpbt27Vxx9/rHHjxmno0KFq1aqVJOnBBx+UyWTSqFGjtGvXLi1btkwLFixQenq6MZbHH39cubm5mjNnjvbs2aPp06fr008/1bhx46r1afoDaxYCiDULAAAA3qj5rVNvuOEGvfPOO5oyZYpmzJihuLg4zZ8/X8OGDTPaTJo0ScePH9djjz2mo0eP6uabb1Zubq4iIiKMNm+99ZbGjRunO+64Q6GhoRo8eLBefPFF47rFYtHatWuVmpqqhIQENW/eXBkZGW73Yrjxxhu1ZMkSTZ06VU899ZSuu+46rVy5Up07d/bh8/BNiJNvo37ncDhksVi0ZcsWNWnSxKe+jh07psTERJWWlrot0sHFuX4GfG64kJqMEeIRl0I8IpjURoz8/Jr3y2xu6GNfp2SxrCDG/YTKQgBRWQAAAPDGaUkhfugD/kKyEEAkCwAAAN4gWQg2LHAGAAAA4BGVhQCisgAAAOANKgvBhmQhgEgWAAAAvEGyEGyYhgQAAADAIyoLAURlAQAAwBsV8r2y4N19FnBxJAsBRLIAAADgDX9MIWIakj8xDQk1IisrS23btlVERIQSExO1devWi7afP3++2rdvr8jISMXGxiotLU0nT56sodGivsvKylKXLl0kSf369SMeUauIRwDBjGQhgFyVBV+Pum7ZsmVKT0/XtGnTtG3bNnXr1k02m00HDx702H7JkiV68sknNW3aNO3evVtvvPGGli1bpqeeeqqGR476yBWPkydPliR17tyZeEStIR6Bc5320wF/IVkIIJKFM+bOnavRo0dr5MiRio+PV3Z2tho1aqRFixZ5bL9p0ybddNNNevDBB9W2bVv1799fDzzwwCV/2wZUhSseH3roIUlnfktLPKK2EI/AuUgWgg3JAgKqvLxchYWFSkpKMs6FhoYqKSlJBQUFHp9z4403qrCw0PjH7+uvv9Z7772nu+6664KvU1ZWJofD4XYA5yIeEUyIRwB1AQucA4gFztLhw4dVUVGh6Ohot/PR0dHas2ePx+c8+OCDOnz4sG6++WY5nU6dPn1aY8aMuWiZPTMzU88884xfx476h3hEMCEeAU/8sZMRuyH5E5WFAGIaUvWsX79eM2fO1Msvv6xt27bp7bff1urVq/Xss89e8DlTpkxRaWmpcezfv78GR4z6jHhEMCEeUf8xDSnYUFlAQDVv3lxhYWEqKSlxO19SUqKYmBiPz/n973+v4cOH69FHH5UkdenSRcePH9djjz2mp59+WqGh5+e44eHhCg8P9/8bQL1ydjx26tTJOE88ojYQjwDqAioLAURlQTKZTEpISFB+fr5xrrKyUvn5+bJarR6f89NPP533D15YWJikuj8tC7WLeEQwIR4BT6gsBBsqCwHEmoUz0tPTlZKSol69eql3796aP3++jh8/rpEjR0qSRowYoV/84hfKzMyUJN19992aO3euevToocTERH355Zf6/e9/r7vvvtv4RxGoLlc8un6Tm5aWRjyi1hCPwLlOS/L1uw9rFvyJZAEBN2TIEB06dEgZGRmy2+3q3r27cnNzjUV9xcXFbr8pmzp1qkJCQjR16lT95z//UYsWLXT33XfrD3/4Q229BdQjrnicOXOmJGnHjh3EI2oN8Qgg2IU468OvroOMw+GQxWJRfn6+Gjdu7FNfx48f1x133KHS0lKZzWY/jbD+c/0M+NxwITUZI8QjLoV4RDCpjRj5+TWvldnsW5XM4aiQxfIVMe4nVBYCiGlIAAAA3qiQ79OQKv0xEPwfFjgDAAAA8IjKQgBRWQAAAPAGlYVgQ2UBAAAAgEdUFgKIygIAAIA3Tsv332VTWfAnkoUAIlkAAADwBslCsGEaEgAAAACPqCwEEJUFAAAAb1BZCDYkCwFEsgAAAOCNCvn+ZZ/vTv5Up6Yhbdy4UXfffbdatWqlkJAQrVy50u36ww8/rJCQELfjzjvvdGtz5MgRDRs2TGazWU2bNtWoUaN07NgxtzafffaZbrnlFkVERCg2NlazZs0K9FsDAAAAgk6dShaOHz+ubt26KSsr64Jt7rzzTn333XfG8de//tXt+rBhw7Rr1y7l5eVp1apV2rhxox577DHjusPhUP/+/dWmTRsVFhbqhRde0PTp0/Xqq696PV5XZcHXAwAA4PJw2k8H/KVOTUMaMGCABgwYcNE24eHhiomJ8Xht9+7dys3N1SeffKJevXpJkl566SXdddddmj17tlq1aqW33npL5eXlWrRokUwmkzp16qSioiLNnTvXLamoCqYhAQAAeOO0pBAf++C7kz/VqcpCVaxfv15RUVFq3769xo4dq++//964VlBQoKZNmxqJgiQlJSUpNDRUW7ZsMdrceuutMplMRhubzaa9e/fqhx9+qLk3AgAAANSyOlVZuJQ777xTgwYNUlxcnL766is99dRTGjBggAoKChQWFia73a6oqCi35zRo0EDNmjWT3W6XJNntdsXFxbm1iY6ONq5deeWV571uWVmZysrKjMcOh0MSlQUAAADvUFkINvUqWRg6dKjx5y5duqhr16669tprtX79et1xxx0Be93MzEw988wz550nWQAAAPCCs9L37/p8dfKrejcN6WzXXHONmjdvri+//FKSFBMTo4MHD7q1OX36tI4cOWKsc4iJiVFJSYlbG9fjC62FmDJlikpLS41j//79/n4rAAAAQI2r18nCv//9b33//fdq2bKlJMlqtero0aMqLCw02qxbt06VlZVKTEw02mzcuFGnTp0y2uTl5al9+/YepyBJZxZVm81mt0NiNyQAAACvVPrpgN/UqWTh2LFjKioqUlFRkSRp3759KioqUnFxsY4dO6aJEydq8+bN+uabb5Sfn6977rlH7dq1k81mkyR17NhRd955p0aPHq2tW7fq448/1rhx4zR06FC1atVKkvTggw/KZDJp1KhR2rVrl5YtW6YFCxYoPT3d6/GSLAAAAHihwk8H/KZOJQuffvqpevTooR49ekiS0tPT1aNHD2VkZCgsLEyfffaZfv3rX+v666/XqFGjlJCQoA8//FDh4eFGH2+99ZY6dOigO+64Q3fddZduvvlmt3soWCwWrV27Vvv27VNCQoJ+97vfKSMjw+ttUwEAAIC6rk4tcL799tsv+pv2NWvWXLKPZs2aacmSJRdt07VrV3344Ydej+9cLHAGAADwgj8qA1QW/KpOJQt1DckCAACAF/yx5oA1C35Vp6YhAQAAAKg5VBYCiMoCAACAF5iGFHRIFgKIZAEAAMALTEMKOkxDAgAAAOARyUIAcZ8FAAAAL1TK93sseFlZmD59ukJCQtyODh06GNdPnjyp1NRUXXXVVWrSpIkGDx6skpIStz6Ki4uVnJysRo0aKSoqShMnTtTp06fd2qxfv149e/ZUeHi42rVrp5ycnPPGkpWVpbZt2yoiIkKJiYnaunWrd28mAEgWAoxEAQAAoIpq6aZsnTp10nfffWccH330kXEtLS1N7777rlasWKENGzbowIEDGjRo0M9DrqhQcnKyysvLtWnTJr355pvKyclRRkaG0Wbfvn1KTk5W3759VVRUpAkTJujRRx912/Z/2bJlSk9P17Rp07Rt2zZ169ZNNptNBw8e9P4N+RHJAgAAAC5rDRo0UExMjHE0b95cklRaWqo33nhDc+fOVb9+/ZSQkKDFixdr06ZN2rx5syRp7dq1+vzzz/WXv/xF3bt314ABA/Tss88qKytL5eXlkqTs7GzFxcVpzpw56tixo8aNG6f77rtP8+bNM8Ywd+5cjR49WiNHjlR8fLyys7PVqFEjLVq0qOY/kLOQLAQQ05AAAAC8UOmnQ5LD4XA7ysrKLviyX3zxhVq1aqVrrrlGw4YNU3FxsSSpsLBQp06dUlJSktG2Q4cOat26tQoKCiRJBQUF6tKli6Kjo402NptNDodDu3btMtqc3YerjauP8vJyFRYWurUJDQ1VUlKS0aa2kCwEEMkCAACAF/w4DSk2NlYWi8U4MjMzPb5kYmKicnJylJubq1deeUX79u3TLbfcoh9//FF2u10mk0lNmzZ1e050dLTsdrskyW63uyUKruuuaxdr43A4dOLECR0+fFgVFRUe27j6qC1snQoAAIB6Z//+/TKbzcbj8PBwj+0GDBhg/Llr165KTExUmzZttHz5ckVGRgZ8nMGOykIAUVkAAADwgh8rC2az2e24ULJwrqZNm+r666/Xl19+qZiYGJWXl+vo0aNubUpKShQTEyNJiomJOW93JNfjS7Uxm82KjIxU8+bNFRYW5rGNq4/aQrIQQCQLAAAAXvDjmoXqOnbsmL766iu1bNlSCQkJatiwofLz843re/fuVXFxsaxWqyTJarVqx44dbrsW5eXlyWw2Kz4+3mhzdh+uNq4+TCaTEhIS3NpUVlYqPz/faFNbSBYAAABw2XriiSe0YcMGffPNN9q0aZPuvfdehYWF6YEHHpDFYtGoUaOUnp6uDz74QIWFhRo5cqSsVqv69OkjSerfv7/i4+M1fPhw/e///q/WrFmjqVOnKjU11ahmjBkzRl9//bUmTZqkPXv26OWXX9by5cuVlpZmjCM9PV2vvfaa3nzzTe3evVtjx47V8ePHNXLkyFr5XFxYsxBA/qgMUFkAAACXjWreJ+G8Przw73//Ww888IC+//57tWjRQjfffLM2b96sFi1aSJLmzZun0NBQDR48WGVlZbLZbHr55ZeN54eFhWnVqlUaO3asrFarGjdurJSUFM2YMcNoExcXp9WrVystLU0LFizQ1Vdfrddff102m81oM2TIEB06dEgZGRmy2+3q3r27cnNzz1v0XNNCnHwb9TuHwyGLxaIlS5aoUaNGPvX1008/6cEHH1RpaanbIh1cnOtnwOeGC6nJGCEecSnEI4JJbcSI8ZrbJHMTH/s6Jll6ihj3E6YhAQAAAPCIaUgBxDQkAAAAL9TCNCRcHMlCAJEsAAAAeIFkIegwDQkAAACAR1QWAojKAgAAgBf8cJ8En58PNyQLAUSyAAAA4AWmIQUdpiEBAAAA8IjKQgBRWQAAAPAClYWgQ7IQQCQLAAAAXmDNQtBhGhIAAAAAj6gsBBCVBQAAAC9UyvdpRFQW/IpkIYBIFgAAALzANKSgwzQkAAAAAB5RWQggKgsAAABeYDekoEOyEEAkCwAAAF4gWQg6TEMCAAAA4BGVhQCisgAAAOAFFjgHHZKFACJZAAAA8ALTkIIO05AAAAAAeERlIcCoDAAAAFQRlYWgQ7IQQExDAgAA8IJTvq854KuTXzENCQAAAIBHVBYCiMoCAACAF5iGFHRIFgKIZAEAAMALbJ0adJiGBAAAAMAjKgsBRGUBAADAC0xDCjokCwFEsgAAAOAFkoWgwzQkAAAAAB5RWQggKgsAAABeYIFz0CFZCCCSBQAAAC8wDSnoMA0JNSIrK0tt27ZVRESEEhMTtXXr1ou2P3r0qFJTU9WyZUuFh4fr+uuv13vvvVdDo0V9l5WVpS5dukiS+vXrRzyiVhGPAIJZnUoWNm7cqLvvvlutWrVSSEiIVq5c6Xbd6XQqIyNDLVu2VGRkpJKSkvTFF1+4tTly5IiGDRsms9mspk2batSoUTp27Jhbm88++0y33HKLIiIiFBsbq1mzZlVrvK7Kgq9HXbds2TKlp6dr2rRp2rZtm7p16yabzaaDBw96bF9eXq5f/vKX+uabb/S3v/1Ne/fu1WuvvaZf/OIXNTxy1EeueJw8ebIkqXPnzsQjag3xCJyjUj9XF6p7MA3Jr+pUsnD8+HF169ZNWVlZHq/PmjVLL774orKzs7VlyxY1btxYNptNJ0+eNNoMGzZMu3btUl5enlatWqWNGzfqscceM647HA71799fbdq0UWFhoV544QVNnz5dr776qtfjJVk4Y+7cuRo9erRGjhyp+Ph4ZWdnq1GjRlq0aJHH9osWLdKRI0e0cuVK3XTTTWrbtq1uu+02devWrYZHjvrIFY8PPfSQJGn+/PnEI2oN8Qico9JPB/ymTiULAwYM0HPPPad77733vGtOp1Pz58/X1KlTdc8996hr167685//rAMHDhgViN27dys3N1evv/66EhMTdfPNN+ull17S0qVLdeDAAUnSW2+9pfLyci1atEidOnXS0KFD9dvf/lZz586tybdab5SXl6uwsFBJSUnGudDQUCUlJamgoMDjc/7xj3/IarUqNTVV0dHR6ty5s2bOnKmKigtPQiwrK5PD4XA7gHMRjwgmxCOAuqBOJQsXs2/fPtntdre/dC0WixITE42/dAsKCtS0aVP16tXLaJOUlKTQ0FBt2bLFaHPrrbfKZDIZbWw2m/bu3asffvjB42tf6C9iKgvS4cOHVVFRoejoaLfz0dHRstvtHp/z9ddf629/+5sqKir03nvv6fe//73mzJmj55577oKvk5mZKYvFYhyxsbF+fR+oH4hHBBPiEfDA1ylI/lggDTf1Jllw/cV6sb907Xa7oqKi3K43aNBAzZo1c2vjqY+zX+NcF/qLmGSheiorKxUVFaVXX31VCQkJGjJkiJ5++mllZ2df8DlTpkxRaWmpcezfv78GR4z6jHhEMCEeUe8xDSnosHWqH0yZMkXp6enGY4fDwW9u/k/z5s0VFhamkpISt/MlJSWKiYnx+JyWLVuqYcOGCgsLM8517NhRdrtd5eXlblUfl/DwcIWHh/t38Kh3zo7HTp06GeeJR9QG4hFAXVBvKguuv1gv9qU0JibmvB0mTp8+rSNHjri18dTH2a9xrvDwcJnNZrdDorIgSSaTSQkJCcrPzzfOVVZWKj8/X1ar1eNzbrrpJn355ZeqrPz5VwP/+te/1LJlS4//EAJVRTwimBCPgAdMQwo69SZZiIuLU0xMjNtfug6HQ1u2bDH+0rVarTp69KgKCwuNNuvWrVNlZaUSExONNhs3btSpU6eMNnl5eWrfvr2uvPJKr8ZEsnBGenq6XnvtNb355pvavXu3xo4dq+PHj2vkyJGSpBEjRmjKlClG+7Fjx+rIkSN6/PHH9a9//UurV6/WzJkzlZqaWltvAfWIKx6XLFkiSUpLSyMeUWuIR+AcJAtBp05NQzp27Ji+/PJL4/G+fftUVFSkZs2aqXXr1powYYKee+45XXfddYqLi9Pvf/97tWrVSgMHDpR0plR75513avTo0crOztapU6c0btw4DR06VK1atZIkPfjgg3rmmWc0atQoTZ48WTt37tSCBQs0b9682njL9cKQIUN06NAhZWRkyG63q3v37srNzTXWghQXFys09Oe8NTY2VmvWrFFaWpq6du2qX/ziF3r88ceNfcgBX7jicebMmZKkHTt2EI+oNcQjgGAX4qxDv7pev369+vbte975lJQU5eTkyOl0atq0aXr11Vd19OhR3XzzzXr55Zd1/fXXG22PHDmicePG6d1331VoaKgGDx6sF198UU2aNDHafPbZZ0pNTdUnn3yi5s2ba/z48V79RexwOGSxWDR79mxFRkb69J5PnDihJ554QqWlpcb0Jlya62fA54YLqckYIR5xKcQjgkltxIjxmrMks29fneQ4IVkmiRj3kzpVWbj99tsvOi0nJCREM2bM0IwZMy7YplmzZka590K6du2qDz/8sNrjdPHHNKI6lMsBAAD4xnUHZ1/7gN/UmzULAAAAAPyLZCGAWOAMAADghVq+z8Lzzz+vkJAQTZgwwTh38uRJpaam6qqrrlKTJk00ePDg83bOLC4uVnJysho1aqSoqChNnDhRp0+fdmuzfv169ezZU+Hh4WrXrp1ycnLOe/2srCy1bdtWERERSkxM1NatW6v/ZvyEZCGASBYAAAC8UIu7IX3yySf605/+pK5du7qdT0tL07vvvqsVK1Zow4YNOnDggAYNGvTzkCsqlJycrPLycm3atElvvvmmcnJylJGRYbTZt2+fkpOT1bdvXxUVFWnChAl69NFHtWbNGqPNsmXLlJ6ermnTpmnbtm3q1q2bbDbbedv+1zSSBQAAAFzWjh07pmHDhum1115z2yq/tLRUb7zxhubOnat+/fopISFBixcv1qZNm7R582ZJ0tq1a/X555/rL3/5i7p3764BAwbo2WefVVZWlsrLyyVJ2dnZiouL05w5c9SxY0eNGzdO9913n9tum3PnztXo0aM1cuRIxcfHKzs7W40aNdKiRYtq9sM4B8lCgFFVAAAAqCI/VhYcDofbUVZWdsGXTU1NVXJyspKSktzOFxYW6tSpU27nO3TooNatW6ugoECSVFBQoC5duhhbHkuSzWaTw+HQrl27jDbn9m2z2Yw+ysvLVVhY6NYmNDRUSUlJRpvaQrIQQExDAgAA8IIf1yzExsbKYrEYR2ZmpseXXLp0qbZt2+bxut1ul8lkUtOmTd3OR0dHy263G23OThRc113XLtbG4XDoxIkTOnz4sCoqKjy2cfVRW+rU1qkAAABAVezfv9/tPgvh4eEe2zz++OPKy8tTRERETQ6vzqCyEEBUFgAAALzgx2lIZrPZ7fCULBQWFurgwYPq2bOnGjRooAYNGmjDhg168cUX1aBBA0VHR6u8vFxHjx51e15JSYliYmIkSTExMeftjuR6fKk2ZrNZkZGRat68ucLCwjy2cfVRW0gWAohkAQAAwAs1vBvSHXfcoR07dqioqMg4evXqpWHDhhl/btiwofLz843n7N27V8XFxbJarZIkq9WqHTt2uO1alJeXJ7PZrPj4eKPN2X242rj6MJlMSkhIcGtTWVmp/Px8o01tYRoSAAAALktXXHGFOnfu7HaucePGuuqqq4zzo0aNUnp6upo1ayaz2azx48fLarWqT58+kqT+/fsrPj5ew4cP16xZs2S32zV16lSlpqYa1YwxY8Zo4cKFmjRpkh555BGtW7dOy5cv1+rVq43XTU9PV0pKinr16qXevXtr/vz5On78uEaOHFlDn4ZnJAsB5I/KAJUFAABw2XDKp5uqGX340bx58xQaGqrBgwerrKxMNptNL7/8snE9LCxMq1at0tixY2W1WtW4cWOlpKRoxowZRpu4uDitXr1aaWlpWrBgga6++mq9/vrrstlsRpshQ4bo0KFDysjIkN1uV/fu3ZWbm3veoueaRrIQQCQLAAAAXvDhpmpuffhg/fr1bo8jIiKUlZWlrKysCz6nTZs2eu+99y7a7+23367t27dftM24ceM0bty4Ko+1JpAsAAAAAPVMeXm5Dh48qMpK91JN69atveqn2gucP/zwQz300EOyWq36z3/+I0n6f//v/+mjjz6qbpf1DgucAQAAvODH+yxcrr744gvdcsstioyMVJs2bRQXF6e4uDi1bdtWcXFxXvdXrcrC//zP/2j48OEaNmyYtm/fbtwRr7S0VDNnzrxkGeZywTQkAAAALwTBNKS67uGHH1aDBg20atUqtWzZUiEhIT71V61k4bnnnlN2drZGjBihpUuXGudvuukmPffccz4NCAAAAED1FBUVqbCwUB06dPBLf9VKFvbu3atbb731vPMWi+W8m1ZczqgsAAAAeIHKgs/i4+N1+PBhv/VXrTULMTEx+vLLL887/9FHH+maa67xeVD1BWsWAAAAvMCaBZ/98Y9/1KRJk7R+/Xp9//33cjgcboe3qlVZGD16tB5//HEtWrRIISEhOnDggAoKCvTEE0/o97//fXW6BAAAAOCjpKQkSWfuTn02p9OpkJAQVVR4V3qpVrLw5JNPqrKyUnfccYd++ukn3XrrrQoPD9cTTzyh8ePHV6fLeolpSAAAAF5gGpLPPvjgA7/2V61kISQkRE8//bQmTpyoL7/8UseOHVN8fLyaNGni18HVdSQLAAAAXqiU71/2L/NpSLfddptf+/Pppmwmk0nx8fH+GgsAAAAAHx09elRvvPGGdu/eLUnq1KmTHnnkEVksFq/7qlaycO+993rcszUkJEQRERFq166dHnzwQbVv37463dcbVBYAAAC84I8Fypd5ZeHTTz+VzWZTZGSkevfuLUmaO3eu/vCHP2jt2rXq2bOnV/1Vazcki8WidevWadu2bQoJCVFISIi2b9+udevW6fTp01q2bJm6deumjz/+uDrd1xvshgQAAOCFCj8dl7G0tDT9+te/1jfffKO3335bb7/9tvbt26df/epXmjBhgtf9VauyEBMTowcffFALFy5UaOiZfKOyslKPP/64rrjiCi1dulRjxozR5MmT9dFHH1XnJQAAAAB46dNPP9Vrr72mBg1+/prfoEEDTZo0Sb169fK6v2pVFt544w1NmDDBSBQkKTQ0VOPHj9err76qkJAQjRs3Tjt37qxO9/UGlQUAAAAvcJ8Fn5nNZhUXF593fv/+/briiiu87q9aycLp06e1Z8+e887v2bPH2Ls1IiLC47qGywnJAgAAgBeYhuSzIUOGaNSoUVq2bJn279+v/fv3a+nSpXr00Uf1wAMPeN1ftaYhDR8+XKNGjdJTTz2lG264QZL0ySefaObMmRoxYoQkacOGDerUqVN1ugcAAABQDbNnz1ZISIhGjBih06dPS5IaNmyosWPH6vnnn/e6v2olC/PmzVN0dLRmzZqlkpISSVJ0dLTS0tI0efJkSVL//v115513Vqf7eoPdkAAAALzATdl8ZjKZtGDBAmVmZuqrr76SJF177bVq1KhRtfqrVrIQFhamp59+Wk8//bQcDoekM/Ojzta6detqDag+IVkAAADwAlun+k2jRo3UpUsXn/vx6aZs0vlJAgAAAICaM2jQIOXk5MhsNmvQoEEXbfv222971Xe1k4W//e1vWr58uYqLi1VeXu52bdu2bdXttt6hMgAAAFBFlfJ9GtFlWFmwWCzGxkLVuUvzxVQrWXjxxRf19NNP6+GHH9bf//53jRw5Ul999ZU++eQTpaam+nWAdRnTkAAAALxQoWru1XlOH5eZxYsXe/yzP1Trx/Hyyy/r1Vdf1UsvvSSTyaRJkyYpLy9Pv/3tb1VaWurXAQIAAAComhMnTuinn34yHn/77beaP3++1q5dW63+qpUsFBcX68Ybb5QkRUZG6scff5R0ZkvVv/71r9UaSH3EfRYAAAC8wE3ZfHbPPffoz3/+syTp6NGj6t27t+bMmaN77rlHr7zyitf9VStZiImJ0ZEjRySd2fVo8+bNkqR9+/bx5fYsJAsAAABe4KZsPtu2bZtuueUWSWfWGMfExOjbb7/Vn//8Z7344ote91etZKFfv376xz/+IUkaOXKk0tLS9Mtf/lJDhgzRvffeW50uAQAAAPjop59+0hVXXCFJWrt2rQYNGqTQ0FD16dNH3377rdf9VWuB86uvvqrKyjM1ntTUVF111VXatGmTfv3rX+s3v/lNdbqsl1jgDAAA4AXus+Czdu3aaeXKlbr33nu1Zs0apaWlSZIOHjxYrVseVCtZCA0NVWjoz0WJoUOHaujQodXpql4jWQAAAPACuyH5LCMjQw8++KDS0tJ0xx13yGq1SjpTZejRo4fX/VX7PgsnT57UZ599poMHDxpVBpdf//rX1e0WAAAAQDXdd999uvnmm/Xdd9+pW7duxvk77rijWssFqpUs5ObmasSIETp8+PB510JCQlRRcZmndP+HygIAAIAXqCz4RUxMjGJiYtzO9e7du1p9VStZGD9+vO6//35lZGQoOjq6Wi98OSBZAAAA8IJTvq85uMy/OvXt29e4m7Mn69at86q/aiULJSUlSk9PJ1EAAAAAgkj37t3dHp86dUpFRUXauXOnUlJSvO6vWsnCfffdp/Xr1+vaa6+tztMvG1QWAAAAvFAh6cK/FK96H5exefPmeTw/ffp0HTt2zOv+qpUsLFy4UPfff78+/PBDdenSRQ0bNnS7/tvf/rY63dY7JAsAAABeIFkImIceeki9e/fW7NmzvXpetZKFv/71r1q7dq0iIiK0fv16t3lRISEhJAsAAABAECkoKFBERITXz6vWevOnn35azzzzjEpLS/XNN99o3759xvH1119Xp0u/mD59ukJCQtyODh06GNdPnjxp3ESuSZMmGjx4sEpKStz6KC4uVnJysho1aqSoqChNnDhRp0+frtZ4XJUFXw8AAIDLQqWfjsvYoEGD3I57771Xffr00ciRI6t18+RqVRbKy8s1ZMgQtxuzBYtOnTrp/fffNx43aPDzW0xLS9Pq1au1YsUKWSwWjRs3ToMGDdLHH38sSaqoqFBycrJiYmK0adMmfffddxoxYoQaNmyomTNnej0WpiEBAAB4gWlI1fb111+rbdu2slgsbudDQ0PVvn17zZgxQ/379/e632olCykpKVq2bJmeeuqp6jw9oBo0aHDevrKSVFpaqjfeeENLlixRv379JEmLFy9Wx44dtXnzZvXp00dr167V559/rvfff1/R0dHq3r27nn32WU2ePFnTp0+XyWSq6bcDAAAAXNJ1112n7777TosXL5YkDRkyRC+++KLPu5dWK1moqKjQrFmztGbNGnXt2vW8Bc5z5871aVC++OKLL9SqVStFRETIarUqMzNTrVu3VmFhoU6dOqWkpCSjbYcOHdS6dWsVFBSoT58+KigoUJcuXdw+VJvNprFjx2rXrl1e3yKbygIAAIAX/DGN6DKdhnTud8Z//vOfOn78uM/9VitZ2LFjh/HFeefOnT4Pwl8SExOVk5Oj9u3b67vvvtMzzzyjW265RTt37pTdbpfJZFLTpk3dnhMdHS273S5Jstvt52VfrseuNp6UlZWprKzMeOxwOCSRLAAAAHiFaUh+46/vkNVKFj744AO/vLi/DRgwwPhz165dlZiYqDZt2mj58uWKjIwM2OtmZmbqmWeeCVj/AAAAwMW4Nvc595yvvEoWBg0adMk2ISEh+p//+Z9qD8ifmjZtquuvv15ffvmlfvnLX6q8vFxHjx51qy6UlJQYaxxiYmK0detWtz5cuyV5WgfhMmXKFKWnpxuPHQ6HYmNjqSwAAAB4o1K+VwYu42lIDz/8sMLDwyWd2QV0zJgxaty4sVu7t99+26t+vUoWzl1dHeyOHTumr776SsOHD1dCQoIaNmyo/Px8DR48WJK0d+9eFRcXy2q1SpKsVqv+8Ic/6ODBg4qKipIk5eXlyWw2Kz4+/oKvEx4ebvxgzkayAAAA4IVK+T4N6TJNFlJSUtweP/TQQ37p16tkwbW6Olg98cQTuvvuu9WmTRsdOHBA06ZNU1hYmB544AFZLBaNGjVK6enpatasmcxms8aPHy+r1ao+ffpIkvr376/4+HgNHz5cs2bNkt1u19SpU5WamuoxGQAAAACCQaC+p1drzUKw+ve//60HHnhA33//vVq0aKGbb75ZmzdvVosWLSRJ8+bNU2hoqAYPHqyysjLZbDa9/PLLxvPDwsK0atUqjR07VlarVY0bN1ZKSopmzJhRrfFQWQAAAPCCPxYns8DZr4Lvrmo+WLp0qQ4cOKCysjL9+9//1tKlS3Xttdca1yMiIpSVlaUjR47o+PHjevvtt89bi9CmTRu99957+umnn3To0CHNnj3b7cZu3uAOzgAAAF6o8NPhhVdeeUVdu3aV2WyW2WyW1WrVP//5T+P6yZMnlZqaqquuukpNmjTR4MGDjTWtLsXFxUpOTlajRo0UFRWliRMn6vTp025t1q9fr549eyo8PFzt2rVTTk7OeWPJyspS27ZtFRERocTExPPW0taGepUsAAAAAN64+uqr9fzzz6uwsFCffvqp+vXrp3vuuUe7du2SJKWlpendd9/VihUrtGHDBh04cMBt05+KigolJyervLxcmzZt0ptvvqmcnBxlZGQYbfbt26fk5GT17dtXRUVFmjBhgh599FGtWbPGaLNs2TKlp6dr2rRp2rZtm7p16yabzaaDBw/W3IfhQYiTX137ncPhkMVi0W9+8xuf7/pcXl6uP/3pTyotLZXZbPbTCOs/18+Azw0XUpMxQjziUohHBJPaiBHjNdtL5jAf+6qQLHvl0/ibNWumF154Qffdd59atGihJUuW6L777pMk7dmzRx07djRu6vvPf/5Tv/rVr3TgwAHj/lzZ2dmaPHmyDh06JJPJpMmTJ2v16tVu9ycbOnSojh49qtzcXEln7hd2ww03aOHChZKkyspKxcbGavz48XryySd9+Uh8QmUhgJiGBAAA4AU/TkNyOBxux9k30L3gy1dUaOnSpTp+/LisVqsKCwt16tQpJSUlGW06dOig1q1bq6CgQJJUUFCgLl26uN3Y12azyeFwGNWJgoICtz5cbVx9lJeXq7Cw0K1NaGiokpKSjDa1hWQBAAAA9U5sbKwsFotxZGZmXrDtjh071KRJE4WHh2vMmDF65513FB8fL7vdLpPJ5HaPLkmKjo6W3W6XJNntdrdEwXXdde1ibRwOh06cOKHDhw+roqLCYxtXH7WlXu2GFGzYDQkAAMALfrzPwv79+92mIV1sG/z27durqKhIpaWl+tvf/qaUlBRt2LDBx4HUDyQLAUSyAAAA4AV/3FDt//pw7W5UFSaTSe3atZMkJSQk6JNPPtGCBQs0ZMgQlZeX6+jRo27VhZKSEmNHzZiYmPN2LXLtlnR2m3N3UCopKZHZbFZkZKTCwsIUFhbmsc25O3fWNKYhAQAAAGeprKxUWVmZEhIS1LBhQ+Xn5xvX9u7dq+LiYlmtVkmS1WrVjh073HYtysvLk9lsVnx8vNHm7D5cbVx9mEwmJSQkuLWprKxUfn6+0aa2UFkIICoLAAAAXqiQ5OtXHy+rE1OmTNGAAQPUunVr/fjjj1qyZInWr1+vNWvWyGKxaNSoUUpPT1ezZs1kNps1fvx4Wa1W9enTR5LUv39/xcfHa/jw4Zo1a5bsdrumTp2q1NRUY+rTmDFjtHDhQk2aNEmPPPKI1q1bp+XLl2v16tXGONLT05WSkqJevXqpd+/emj9/vo4fP66RI0f6+IH4hspCALEb0s+qe5ORpUuXKiQkRAMHDgzsAHFZycrKUpcuXSRJ/fr1Ix5Rq4hH4CyVfjq8cPDgQY0YMULt27fXHXfcoU8++URr1qzRL3/5S0nSvHnz9Ktf/UqDBw/WrbfeqpiYGL399tvG88PCwrRq1SqFhYXJarXqoYce0ogRIzRjxgyjTVxcnFavXq28vDx169ZNc+bM0euvvy6bzWa0GTJkiGbPnq2MjAx1795dRUVFys3NPW/Rc03jPgsB4Nor+JFHHvHLfRYWLVpUp/fDXrZsmUaMGKHs7GwlJiZq/vz5WrFihfbu3auoqKgLPu+bb77RzTffrGuuuUbNmjXTypUrq/ya7COOC3HF47x585SamqqUlBT9/e9/Jx5RK4hHBKNavc/CLySzj7/KdlRKlv/4dp8F/IzKQgBRWThj7ty5Gj16tEaOHKn4+HhlZ2erUaNGWrRo0QWfU1FRoWHDhumZZ57RNddcU4OjRX3niseHHnpIkjR//nziEbWGeATO4cf7LMA/SBYCiGSh+jcZmTFjhqKiojRq1KgqvU5ZWdl5N18BzkU8IpgQj4AHJAtBh2QBAVWdm4x89NFHeuONN/Taa69V+XUyMzPdbrwSGxvr07hRPxGPCCbEI4C6gGQhgKgseO/HH3/U8OHD9dprr6l58+ZVft6UKVNUWlpqHPv37w/gKHG5IB4RTIhHXBZqYYEzLo6tUwOIrVOl5s2be3WTka+++krffPON7r77buNcZeWZ/+sbNGigvXv36tprrz3veeHh4Re9MyMgucdjp06djPPEI2oD8Qh4UCnft06t21+dgg6VBQSUtzcZ6dChg3bs2KGioiLj+PWvf62+ffuqqKiI8jl8QjwimBCPAOoCKgsBRGXhjEvdZGTEiBH6xS9+oczMTEVERKhz585uz3fdXv3c80B1uOLR9ZvctLQ04hG1hngEzlEpKcTHPur+V6egQrIQQCQLZwwZMkSHDh1SRkaG7Ha7unfv7naTkeLiYoWGUuRCzXDF48yZMyVJO3bsIB5Ra4hH4BwVIlkIMtyULQBcNxZ56KGH/HJTtr/85S/cWMRL3HQIl1KTMUI84lKIRwSTWr0pWxPJ7GOy4HBKlmPclM1fqCwEEJUFAAAAL1BZCDokCwFEsgAAAOAF1iwEHSZCAgAAAPCIykIAUVkAAADwAtOQgg7JQgCRLAAAAHiBZCHoMA0JAAAAgEdUFgKIygIAAIAXnKIyEGRIFgKIZAEAAKDqKv7v8LUP+A/TkAAAAAB4RGUhwKgMAAAAVA2VheBDshBATEMCAACousr/O3ztA/7DNCQAAAAAHlFZCCAqCwAAAFXHNKTgQ7IQQCQLAAAAVcc0pODDNCQAAAAAHlFZCCAqCwAAAFXHNKTgQ7IQQCQLAAAAVVcp37/sMw3Jv5iGBAAAAMAjKgsBRGUBAACg6ljgHHxIFgKIZAEAAKDqWLMQfJiGBAAAAMAjKgsBRGUBAACg6qgsBB+ShQAiWQAAAKg61iwEH6YhAQAAAPCIykIAUVkAAACoOqYhBR+ShQAiWQAAAKg6piEFH6YhAQAAAPCIykIAUVkAAACoukr5Po2IyoJ/UVm4gKysLLVt21YRERFKTEzU1q1bve7DlSz4egAAAFwOKvx0wH9IFjxYtmyZ0tPTNW3aNG3btk3dunWTzWbTwYMHa3toAAAAQI0hWfBg7ty5Gj16tEaOHKn4+HhlZ2erUaNGWrRokVf9UFkAAACouko/HfAf1iyco7y8XIWFhZoyZYpxLjQ0VElJSSooKPD4nLKyMpWVlRmPHQ6HJNYsAAAAeIOtU4MPlYVzHD58WBUVFYqOjnY7Hx0dLbvd7vE5mZmZslgsxhEbG1sTQwUAAAACimTBD6ZMmaLS0lLj2L9/vySmIQEAAHiDBc7Bh2ThHM2bN1dYWJhKSkrczpeUlCgmJsbjc8LDw2U2m90OiWQBAADAG7WxZiEzM1M33HCDrrjiCkVFRWngwIHau3evW5uTJ08qNTVVV111lZo0aaLBgwef912xuLhYycnJatSokaKiojRx4kSdPn3arc369evVs2dPhYeHq127dsrJyTlvPP7YkdOfSBbOYTKZlJCQoPz8fONcZWWl8vPzZbVaa3FkAAAA8LcNGzYoNTVVmzdvVl5enk6dOqX+/fvr+PHjRpu0tDS9++67WrFihTZs2KADBw5o0KBBxvWKigolJyervLxcmzZt0ptvvqmcnBxlZGQYbfbt26fk5GT17dtXRUVFmjBhgh599FGtWbPGaBOMO3KGOPnV9XmWLVumlJQU/elPf1Lv3r01f/58LV++XHv27DlvLYMnDodDFotFSUlJatiwoU9jOXXqlN5//32VlpYaFQtcmutnwOeGC6nJGCEecSnEI4JJbcSI6zXXSWriY1/HJPWTqj3+Q4cOKSoqShs2bNCtt96q0tJStWjRQkuWLNF9990nSdqzZ486duyogoIC9enTR//85z/1q1/9SgcOHDC+K2ZnZ2vy5Mk6dOiQTCaTJk+erNWrV2vnzp3Gaw0dOlRHjx5Vbm6uJCkxMVE33HCDFi5cKOnML6xjY2M1fvx4Pfnkk759MNVEZcGDIUOGaPbs2crIyFD37t1VVFSk3NzcKiUK52IKEgAAQNU45fsUJNe3J4fD4XacvXPlxZSWlkqSmjVrJkkqLCzUqVOnlJSUZLTp0KGDWrdubeyUWVBQoC5durh9V7TZbHI4HNq1a5fR5uw+XG1cfbh25Dy7zaV25KwJJAsXMG7cOH377bcqKyvTli1blJiYWNtDAgAAQBXFxsa67VaZmZl5yedUVlZqwoQJuummm9S5c2dJkt1ul8lkUtOmTd3anr1Tpt1u97iTpuvaxdo4HA6dOHGiWjty1gTusxBA3GcBAACg6vx5n4X9+/e7TUMKDw+/5HNTU1O1c+dOffTRRz6Oov4gWQAAAEC9c/YOlVUxbtw4rVq1Shs3btTVV19tnI+JiVF5ebmOHj3qVl04e6fMmJiY83Ytcu2WdHYbT7ttms1mRUZGKiwszOsdOWsC05ACiK1TAQAAqq427rPgdDo1btw4vfPOO1q3bp3i4uLcrickJKhhw4ZuO2Xu3btXxcXFxk6ZVqtVO3bscNu1KC8vT2azWfHx8Uabs/twtXH1Eaw7clJZCCCmIQEAAFRdde6T4KkPb6SmpmrJkiX6+9//riuuuMJYH2CxWBQZGSmLxaJRo0YpPT1dzZo1k9ls1vjx42W1WtWnTx9JUv/+/RUfH6/hw4dr1qxZstvtmjp1qlJTU43pT2PGjNHChQs1adIkPfLII1q3bp2WL1+u1atXG2NJT09XSkqKevXqZezIefz4cY0cOdLHT6X6SBYAAABw2XrllVckSbfffrvb+cWLF+vhhx+WJM2bN0+hoaEaPHiwysrKZLPZ9PLLLxttw8LCtGrVKo0dO1ZWq1WNGzdWSkqKZsyYYbSJi4vT6tWrlZaWpgULFujqq6/W66+/LpvNZrQZMmSIDh06pIyMDNntdnXv3r3aO3L6C/dZCADXXsG33367GjTwLR87ffq01q9fz37YXmIfcVwK+9ojmBCPCCa1eZ+Ff0hq7GNfxyX9WtW/zwLcUVkIIKYhAQAAVJ0/d0OCf7DAGQAAAIBHVBYCiMoCAABA1dXGAmdcHMlCAJEsAAAAVF2lfJ9GRLLgX0xDAgAAAOARlYUAorIAAABQdUxDCj4kCwFEsgAAAFB17IYUfJiGBAAAAMAjKgsBRGUBAACg6qgsBB+ShQAiWQAAAKg61iwEH6YhAQAAAPCIykIAUVkAAACoOqYhBR+ShQAiWQAAAKg6koXgwzQkAAAAAB5RWQggKgsAAABV55TvC5T55uRfJAsBRLIAAABQdUxDCj5MQwIAAADgEZWFAKKyAAAAUHXcZyH4kCwEEMkCAABA1TENKfgwDQkAAACAR1QWAozKAAAAQNVQWQg+VBYCyDUNydejPsjKylLbtm0VERGhxMREbd269YJtX3vtNd1yyy268sordeWVVyopKemi7QFvZWVlqUuXLpKkfv36EY+oVcQj8LNKPx3wH5IFBNyyZcuUnp6uadOmadu2berWrZtsNpsOHjzosf369ev1wAMP6IMPPlBBQYFiY2PVv39//ec//6nhkaM+csXj5MmTJUmdO3cmHlFriEcAwS7EWV9+dR1EHA6HLBaLevToobCwMJ/6qqio0Pbt21VaWiqz2eynEdasxMRE3XDDDVq4cKEkqbKyUrGxsRo/fryefPLJSz6/oqJCV155pRYuXKgRI0ZU6TVdP4O6/LkhMFzxOHPmTFksFv3www/q1KkT8YhaQTwiGNVGjLhec5akSB/7OiFpkkSM+wmVhQBiGpJUXl6uwsJCJSUlGedCQ0OVlJSkgoKCKvXx008/6dSpU2rWrNkF25SVlcnhcLgdwLmIRwQT4hE4X6V+XrdQ3YNpSP5FsoCAOnz4sCoqKhQdHe12Pjo6Wna7vUp9TJ48Wa1atXL7B/VcmZmZslgsxhEbG+vTuFE/EY8IJsQjgLqAZCGAqCz47vnnn9fSpUv1zjvvKCIi4oLtpkyZotLSUuPYv39/DY4SlwviEcGEeER9xALn4MPWqQHETdmk5s2bKywsTCUlJW7nS0pKFBMTc9Hnzp49W88//7zef/99de3a9aJtw8PDFR4e7vN4Ub+dHY+dOnUyzhOPqA3EI3A+tk4NPlQWEFAmk0kJCQnKz883zlVWVio/P19Wq/WCz5s1a5aeffZZ5ebmqlevXjUxVFwGiEcEE+IRQF1AZSGAqCyckZ6erpSUFPXq1Uu9e/fW/Pnzdfz4cY0cOVKSNGLECP3iF79QZmamJOmPf/yjMjIytGTJErVt29aYu9ukSRM1adKk1t4H6gdXPLp+k5uWlkY8otYQj4A7f0wjYhqSf5EsBBDJwhlDhgzRoUOHlJGRIbvdru7duys3N9dY1FdcXKzQ0J+LXK+88orKy8t13333ufUzbdo0TZ8+vSaHjnrIFY8zZ86UJO3YsYN4RK0hHgF3TEMKPtxnIQBcewV37tzZL/dZ2LlzJ3sFe4l9xHEpNRkjxCMuhXhEMKnN+yxkSLrwcv2qOSlphrjPgr9QWQggKgsAAABVR2Uh+JAsBBDJAgAAQNWxZiH4sBsSAAAAAI+oLAQQlQUAAICqq5Tv04ioLPgXyUIAkSwAAABUHWsWgg/TkAAAAAB4RGUhgKgsAAAAVB0LnIMPyUIAkSwAAABUHdOQgk+9mobUtm1bhYSEuB3PP/+8W5vPPvtMt9xyiyIiIhQbG6tZs2ad18+KFSvUoUMHRUREqEuXLnrvvfdq6i0AAAAAQaNeJQuSNGPGDH333XfGMX78eOOaw+FQ//791aZNGxUWFuqFF17Q9OnT9eqrrxptNm3apAceeECjRo3S9u3bNXDgQA0cOFA7d+70eiyuyoKvBwAAwOWg0k8H/KfeTUO64oorFBMT4/HaW2+9pfLyci1atEgmk0mdOnVSUVGR5s6dq8cee0yStGDBAt15552aOHGiJOnZZ59VXl6eFi5cqOzsbK/GwjQkAACAqmMaUvCpd5WF559/XldddZV69OihF154QadPnzauFRQU6NZbb5XJZDLO2Ww27d27Vz/88IPRJikpya1Pm82mgoKCmnkDAAAAQJCoV8nCb3/7Wy1dulQffPCBfvOb32jmzJmaNGmScd1utys6OtrtOa7Hdrv9om1c1z0pKyuTw+FwOySmIQEAAHijwk+HNzZu3Ki7775brVq1UkhIiFauXOl23el0KiMjQy1btlRkZKSSkpL0xRdfuLU5cuSIhg0bJrPZrKZNm2rUqFE6duyYW5u6um426JOFJ5988rxFy+cee/bskSSlp6fr9ttvV9euXTVmzBjNmTNHL730ksrKygI6xszMTFksFuOIjY2VRLIAAADgDad8X6/g7Ten48ePq1u3bsrKyvJ4fdasWXrxxReVnZ2tLVu2qHHjxrLZbDp58qTRZtiwYdq1a5fy8vK0atUqbdy40ZjiLtX8ull/Cvo1C7/73e/08MMPX7TNNddc4/F8YmKiTp8+rW+++Ubt27dXTEyMSkpK3Nq4HrvWOVyozYXWQUjSlClTlJ6ebjx2OBxGwgAAAIDgNWDAAA0YMMDjNafTqfnz52vq1Km65557JEl//vOfFR0drZUrV2ro0KHavXu3cnNz9cknn6hXr16SpJdeekl33XWXZs+erVatWtX4ull/CvrKQosWLdShQ4eLHmevQThbUVGRQkNDFRUVJUmyWq3auHGjTp06ZbTJy8tT+/btdeWVVxpt8vPz3frJy8uT1Wq94BjDw8NlNpvdDheqCgAAAFXjz2lI504Rr85Mk3379slut7utZ7VYLEpMTDTWsxYUFKhp06ZGoiBJSUlJCg0N1ZYtW4w2dXXdbNAnC1VVUFCg+fPn63//93/19ddf66233lJaWpoeeughIxF48MEHZTKZNGrUKO3atUvLli3TggUL3KoCjz/+uHJzczVnzhzt2bNH06dP16effqpx48Z5PSamIQEAAFSdP5OF2NhYt2nimZmZXo/HtWb1YutZ7Xa78YtplwYNGqhZs2aXXBN79mtUZ91sTQj6aUhVFR4erqVLl2r69OkqKytTXFyc0tLS3BIBi8WitWvXKjU1VQkJCWrevLkyMjLc5pTdeOONWrJkiaZOnaqnnnpK1113nVauXKnOnTvXxtsCAABANezfv99ttkd4eHgtjqbuqjfJQs+ePbV58+ZLtuvatas+/PDDi7a5//77df/99/s8Ju6zAAAAUHX+uKma6/nnTg2vDtea1ZKSErVs2dI4X1JSou7duxttDh486Pa806dP68iRI5dcE3v2a1Rn3WxNqDfTkIIR05AAAACqrja2Tr2YuLg4xcTEuK1ndTgc2rJli7Ge1Wq16ujRoyosLDTarFu3TpWVlUpMTDTaBGLdbE0gWQAAAMBl69ixYyoqKlJRUZGkM4uai4qKVFxcrJCQEE2YMEHPPfec/vGPf2jHjh0aMWKEWrVqpYEDB0qSOnbsqDvvvFOjR4/W1q1b9fHHH2vcuHEaOnSoWrVqJanm1836U72ZhhSMmIYEAABQdf6chlRVn376qfr27Ws8dn2BT0lJUU5OjiZNmqTjx4/rscce09GjR3XzzTcrNzdXERERxnPeeustjRs3TnfccYdCQ0M1ePBgvfjii8b1urxuNsTJt1G/czgcxs3ZQkN9K95UVlZq//79Ki0t9Xne3eXE9TPgc8OF1GSMEI+4FOIRwaQ2YsT1mv8lyfOG+FVXLmm5RIz7CdOQAAAAAHjENKQAYhoSAABA1VXK9wXKvk5jgjuShQAiWQAAAKi62lizgItjGhIAAAAAj6gsBBCVBQAAgKqrkO+/yfbnfRZAshBQJAsAAABVR7IQfJiGBAAAAMAjKgsBRGUBAACg6ljgHHxIFgKIZAEAAKDqmIYUfJiGBAAAAMAjKgsBRGUBAACg6piGFHxIFgKIZAEAAKDquINz8GEaEgAAAACPqCwEEJUFAACAqquQFOKHPuA/JAsBRLIAAABQdaxZCD5MQwIAAADgEZWFAKKyAAAAUHVMQwo+JAsBRLIAAABQdSQLwYdpSAAAAAA8orIQYFQGAAAAqoYFzsGHZCGA/JEokGwAAIDLBdOQgg/TkAAAAAB4RGUhgKgsAAAAVJ1Tvk8j4puTf5EsBBDJAgAAQNX5YwoR05D8i2lIAAAAADyishBAVBYAAACqjspC8CFZCCCSBQAAgKqrlO+7IbF1qn8xDQkAAACAR1QWAojKAgAAQNUxDSn4kCwEEMkCAABA1ZEsBB+mIQEAAADwiMpCAFFZAAAAqDoWOAcfkoUAIlkAAACoOn980SdZ8C+mIQEAAADwiMpCAFFZAAAAqDoqC8GHZCGASBYAAACqrkKSr998SBb8i2lIAAAAADyishBAVBYAAACqjspC8CFZCCCSBQAAgKpjzULwYRoSAAAAAI+oLAQQlQUAAICqYxpS8CFZCCCSBQAAgKqrlO/JAt+c/ItpSAAAAAA8qjPJwh/+8AfdeOONatSokZo2beqxTXFxsZKTk9WoUSNFRUVp4sSJOn36tFub9evXq2fPngoPD1e7du2Uk5NzXj9ZWVlq27atIiIilJiYqK1bt1ZrzE6n0y9HfeDtZ7pixQp16NBBERER6tKli957770aGikuB1lZWerSpYskqV+/fsQjahXxCPys0k9Hdfjr+199U2eShfLyct1///0aO3asx+sVFRVKTk5WeXm5Nm3apDfffFM5OTnKyMgw2uzbt0/Jycnq27evioqKNGHCBD366KNas2aN0WbZsmVKT0/XtGnTtG3bNnXr1k02m00HDx70eswkC2d4+5lu2rRJDzzwgEaNGqXt27dr4MCBGjhwoHbu3FnDI0d95IrHyZMnS5I6d+5MPKLWEI+Auwo/Hd7y5/e/esdZxyxevNhpsVjOO//ee+85Q0NDnXa73Tj3yiuvOM1ms7OsrMzpdDqdkyZNcnbq1MnteUOGDHHabDbjce/evZ2pqanG44qKCmerVq2cmZmZVR5jaWmpU2emzDlDQkJ8Olz9lJaWVvn1g423n+l//dd/OZOTk93OJSYmOn/zm99U+TVdP4O6/LkhMFzx6IqRH374gXhErSEeEYxqI0Zcr9lEcl7h49GkGt+d/PH9r76qNwucCwoK1KVLF0VHRxvnbDabxo4dq127dqlHjx4qKChQUlKS2/NsNpsmTJgg6Uz1orCwUFOmTDGuh4aGKikpSQUFBRd87bKyMpWVlRmPHQ6H8WdnPagM+KI6n2lBQYHS09PdztlsNq1cufKCr3Puz6C0tFSS+88CcMXj448/bsRGSEgI8YhaQTwiWLlioza+w1RICvGxD9eoz43x8PBwhYeHn9e+ut//Lhf1Jlmw2+1uiYIk47Hdbr9oG4fDoRMnTuiHH35QRUWFxzZ79uy54GtnZmbqmWeecTvXpEkTHTt2rNrv52wxMTEymUx+6aumHT582OvP9EI/J9fP0RNPPwNJio2NrcaoUd8NGzbM+PP3339PPKJWEY8IVt9//70sFkuNvJbJZFJMTMxFY9kbTZo0OS/Gp02bpunTp5/XtjrfVS4ntZosPPnkk/rjH/940Ta7d+9Whw4damhE1TNlypTzftPjdDoVEuJrbnyGyWRSRESEX/qqr879GRw9elRt2rRRcXFxjf1FF8wcDodiY2O1f/9+mc3m2h5Orfnuu+/UoUMH5eXlqX379mrdurWaNWvm99chHi+OeDyDeAwOxOP5SktLAxaPFxIREaF9+/apvLzcL/15+h7mqaqAS6vVZOF3v/udHn744Yu2ueaaa6rUV0xMzHmr1ktKSoxrrv+6zp3dxmw2KzIyUmFhYQoLC/PYxtWHJxcqa0Fq3ry515/phX5O1fkZWCwW/vI/i9lsvqw/j4iICIWFhenYsWPGl6TQ0FDisZYQj8RjMLnc49GT0NCa3QcnIiKiVn45Wp3vKpeTWt0NqUWLFurQocNFj6pOv7FardqxY4fbqvW8vDyZzWbFx8cbbfLz892el5eXJ6vVKunMb/ATEhLc2lRWVio/P99oA+9U5zO91M8JqC7iEcGEeASCA9//LqF211dX3bfffuvcvn2785lnnnE2adLEuX37duf27dudP/74o9PpdDpPnz7t7Ny5s7N///7OoqIiZ25urrNFixbOKVOmGH18/fXXzkaNGjknTpzo3L17tzMrK8sZFhbmzM3NNdosXbrUGR4e7szJyXF+/vnnzscee8zZtGlTt12W4J1LfabDhw93Pvnkk0b7jz/+2NmgQQPn7Nmznbt373ZOmzbN2bBhQ+eOHTuq/Jrs9uGOz+Nnrnh85ZVXnJKcDz/8MPFYw/g8fkY81j4+j/Ndjp8J3/8urM4kCykpKcY2omcfH3zwgdHmm2++cQ4YMMAZGRnpbN68ufN3v/ud89SpU279fPDBB87u3bs7TSaT85prrnEuXrz4vNd66aWXnK1bt3aaTCZn7969nZs3bw7wu6v/LvaZ3nbbbc6UlBS39suXL3def/31TpPJ5OzUqZNz9erVXr3eyZMnndOmTXOePHnSH8Ov8/g83L300kvO2NhYZ1hYmPOGG24gHmsYn4c74rF28Xmc73L9TPj+51mI03mZ7+0JAAAAwKM6cwdnAAAAADWLZAEAAACARyQLAAAAADwiWQAAAADgEckC6qysrCy1bdtWERERSkxMPO+mfOdasWKFOnTooIiICHXp0kXvvfdeDY20ZnjzeeTk5CgkJMTtqE93Cd+4caPuvvtutWrVSiEhIVq5cuUln7N+/Xr17NlT4eHhateunXJycrx6TeLRHfH4M+Kx9hGPP6uNeETdRrKAOmnZsmVKT0/XtGnTtG3bNnXr1k02m83tpnxn27Rpkx544AGNGjVK27dv18CBAzVw4EDt3LmzhkceGN5+HtKZu5V+9913xvHtt9/W4IgD6/jx4+rWrZuysrKq1H7fvn1KTk5W3759VVRUpAkTJujRRx/VmjVrqvR84tEd8eiOeKxdxKO7mo5H1AO1vXcrUB29e/d2pqamGo8rKiqcrVq1cmZmZnps/1//9V/O5ORkt3OJiYnO3/zmNwEdZ03x9vNYvHix02Kx1NDoapck5zvvvHPRNpMmTXJ26tTJ7dyQIUOcNputSq9BPLojHi+MeKx5xOOF1UQ8ou6jsoA6p7y8XIWFhUpKSjLOhYaGKikpSQUFBR6fU1BQ4NZekmw22wXb1yXV+Twk6dixY2rTpo1iY2N1zz33aNeuXTUx3KDkS3wQj+6IR98Rj/5DPPquPscHqoZkAXXO4cOHVVFRoejoaLfz0dHRstvtHp9jt9u9al+XVOfzaN++vRYtWqS///3v+stf/qLKykrdeOON+ve//10TQw46F4oPh8OhEydOXPS5xKM74tF3xKP/EI++8yUeUT80qO0BAKh5VqtVVqvVeHzjjTeqY8eO+tOf/qRnn322FkeGyxHxiGBCPALuqCygzmnevLnCwsJUUlLidr6kpEQxMTEenxMTE+NV+7qkOp/HuRo2bKgePXroyy+/DMQQg96F4sNsNisyMvKizyUe3RGPviMe/Yd49J0v8Yj6gWQBdY7JZFJCQoLy8/ONc5WVlcrPz3f7bdDZrFarW3tJysvLu2D7uqQ6n8e5KioqtGPHDrVs2TJQwwxqvsQH8eiOePQd8eg/xKPv6nN8oIpqe4U1UB1Lly51hoeHO3Nycpyff/6587HHHnM2bdrUabfbnU6n0zl8+HDnk08+abT/+OOPnQ0aNHDOnj3buXv3bue0adOcDRs2dO7YsaO23oJfeft5PPPMM841a9Y4v/rqK2dhYaFz6NChzoiICOeuXbtq6y341Y8//ujcvn27c/v27U5Jzrlz5zq3b9/u/Pbbb51Op9P55JNPOocPH260//rrr52NGjVyTpw40bl7925nVlaWMywszJmbm1ul1yMe3RGP7ojH2kU8uqvpeETdR7KAOuull15ytm7d2mkymZy9e/d2bt682bh22223OVNSUtzaL1++3Hn99dc7TSaTs1OnTs7Vq1fX8IgDy5vPY8KECUbb6Oho51133eXctm1bLYw6MD744AOnpPMO12eQkpLivO222857Tvfu3Z0mk8l5zTXXOBcvXuzVaxKP7ojHnxGPtY94/FltxCPqthCn0+ms6WoGAAAAgODHmgUAAAAAHpEsAAAAAPCIZAEAAACARyQLAAAAADwiWQAAAADgEckCAAAAAI9IFgAAAAB4RLIAAAAAwCOSBQAAAAAekSwAAAAA8IhkAQAAAIBHJAsAAAAAPPr/b0bOVd1KnN0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x500 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c41848b7b584993bf77a9de53b89e7e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=36, description='sx', max=71), IntSlider(value=62, description='sy', max…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# display \n",
    "show_nparray_3d_image(ct_array, image2=dosemap_array, rot90=(2, 2, 0), figsize=(10, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "WARNING : the computed dose values are not really realistic. We consider here that there is no washout, i.e. no biological decay, only the physical decay. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-16T14:46:34.735171Z",
     "start_time": "2024-09-16T14:46:34.722273Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ROI image size = (124, 72, 190)\n",
      "Mean dose in the liver = 22.13662370173907 Gy\n"
     ]
    }
   ],
   "source": [
    "# read some organs segmentation (the masks have the same size as the ct)\n",
    "roi_filename = data_folder / \"liver.nii.gz\"\n",
    "roi = sitk.ReadImage(roi_filename)\n",
    "roi_array = sitk.GetArrayViewFromImage(roi)\n",
    "print('ROI image size =', roi.GetSize())\n",
    "\n",
    "# compute the mean dose in the region\n",
    "mean_dose = np.mean(dosemap_array[roi_array == 1])\n",
    "print('Mean dose in the liver =', mean_dose, 'Gy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "Compute the dose in others ROIs. \n"
   ]
  }
 ],
 "metadata": {
  "deepnote_execution_queue": [],
  "deepnote_notebook_id": "97ba8bac-eb79-4668-ab41-70a851c67caa",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
